##### stratification (correction en fonction de l'origine) #####

# creation d'un tableau de l'ensemble des donnees (avec hardy-weinberg, pourcentage de génotypage et maf)
source("GWA_analysis_functions.R")
create_summary_data(inbfile)

# Suppression des SNPs selon leur équilibre de hardy-weinberg, maff, missing
snp_selection(HW,miss,maf)

system(paste("plink --noweb --bfile ",inbfile," --exclude exclusionSNP.txt --make-bed --out forStrat",sep=""))

# Suppression des SNPs en fonction de leur fréquence chez les cas et les contrôles 
del_freq(inbfile="forStrat",th_freq=0.05)
system(paste("plink --noweb --bfile ",inbfile," --exclude SNPoutfreq.txt --make-bed --out forStrat",sep=""))

# Sélection des SNPs indépendants
system(paste("plink --noweb --bfile ",inbfile," --hwe 1e-04 --geno 0.05 --maf 0.1 --mind 0.05 --indep-pairwise 50 10 0.2",sep=""))
system(paste("plink --noweb --bfile ",inbfile," --extract plink.prune.in --make-bed --out forStrat",sep=""))

# Elimination des SNPs en LD résiduel
system(paste("plink --noweb --bfile forStrat --r2 --ld-window 50 --ld-window-kb 5000 --ld-window-r2 0.2 --out ld",sep=""))
system(paste("perl snpIn.pl ld.ld",threshold))
system(paste("plink --noweb --bfile forStrat --exclude SNP_out.txt --make-bed --out forStrat2",sep=""))
#verification
system(paste("plink --noweb --bfile forStrat2 --r2 --out verif",sep=""))

# Détermination de la matrice IBD
system(paste("nohup plink --noweb --bfile forStrat2 --genome &",sep="")) 

PIHAT.plots()
PIHAT.exclusionlist(genomefile="plink.genome",famfile="forStrat2.fam",th_indiv=0.1,th_mean=0.015,controls=c("ADN103M2","ADN103M3","ADN103M4","ADN103M5","ADN103M6"))

# exclusion des individus ayant un trop fort pi-hat moyen
system(paste("plink --noweb --bfile forStrat2 --remove exclude_indiv_meanpi-hat.txt --make-bed --out forStrat2",sep=""))

# exclusion d'individu par rapport au PIHAT
system(paste("plink --noweb --bfile forStrat2 --remove list-indiv-pihat-sup0.1.txt --make-bed --out forStrat2",sep=""))

# Détermination des composantes principales
system(paste("plink --noweb --bfile forStrat2 --read-genome plink.genome --mds-plot 4 --cluster --out strat",sep=""))

# création d'un fichier de phenotype à partir des covariables
create_strat_phenofile()

# création de plots montrant l'association des individus aux composantes
plot.indiv.components(nb_comp=4)

# creation de la liste d'individus à éliminer (en fonction des associations aux composantes)
# a faire dans R à partir du fichier pheno_strat.txt

# elimination des individus en fonction des associations aux composantes
system((paste("plink --noweb --bfile forStrat2 --remove indiv-components.txt --make-bed --out forStrat2",sep=""))

# Détermination des composantes principales
system(paste("plink --noweb --bfile forStrat2 --read-genome plink.genome --mds-plot 4 --cluster --out strat",sep=""))

# création d'un fichier de phenotype à partir des covariables
create_strat_phenofile()

# analyse d'association avec ces covariables 
system(paste("plink --noweb --bfile ",inbfile," --pheno pheno_strat.txt --all-pheno --assoc --out pheno",sep=""))

# création d'un tableau récapitulatif des statistiques par rapport aux composantes
summary.qassoc()

qqplot_components()
